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The Quantum Monte Carlo method for spin 1/2 fermions at finite temperature is formulated 
for dilute systems with an s-wave interaction. The motivation and the formalism are discussed 
■ along with descriptions of the algorithm and various numerical issues. We report on results for the 

' energy, entropy and chemical potential as a function of temperature. We give upper bounds on the 

, critical temperature Tc for the onset of superfluidity, obtained by studying the finite size scaling of 

the condensate fraction. All of these quantities were computed for couplings around the unitary 
regime in the range -0.5 < {kpay^ < 0.2, where a is the s-wave scattering length and kp is the 
Fermi momentum of a non-interacting gas at the same density. In all cases our data is consistent 
with normal Fermi gas behavior above a characteristic temperature To > Tc, which depends on the 
coupling and is obtained by studying the deviation of the caloric curve from that of a free Fermi 
gas. For Tc < T < To we find deviations from normal Fermi gas behavior that can be attributed to 
pairing effects. Low temperature results for the energy and the pairing gap are shown and compared 
' with Green Function Monte Carlo results by other groups. 
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c5 , I. INTRODUCTION 

<— > 

, The last few years have witnessed an extraordinary progress in the field of cold fermionic atoms, particularly since 

■ the experimental observation of superfluidity in these systems [l|. Ultracold atomic gases provide an exceptional 
G opportunity to investigate strongly correlated fermions. Taking advantage of Feshbach resonances, experimentalists 

vary the strength of the interaction between the atoms at will, offering an unprecedented ability to study the BCS-BEC 
crossover. It is by now well established that, for the case of broad resonances, the physics of these systems is chiefly 
captured by a model of spin 1/2 fermions with a contact s-wave attractive interaction. On the theoretical side, our 
^ ' overall understanding of these remarkable many-body systems has improved dramatically, although many questions 
still remain unanswered. Systems in the BCS-BEC crossover are strongly correlated, and nonperturbative approaches 
are needed. As a consequence, numerical simulations of spin 1/2 fermions at zero and non-zero temperature have 
lately attracted extraordinary attention (see refs.jl, H, 0, 0, H, 0, H, Q)- 

Of particular interest within the BCS-BEC crossover is the so-called unitary regime (see [13, [H [il [il). This 
, regime is formally deflned as the limit of diluteness with respect to the range of the interaction tq , and large scattering 
' length a, such that nr^ <C 1 <C n|a|^, where n is the particle number density. The thermal behavior in the crossover 
[ is characterized by a dimensionless universal function conventionally called ^(T/ej?, l/fci?a), that depends on the 

■ temperature T (in units of the free gas Fermi energy ep — h^kp/2m, where kp = (Stt^tt,)^/'^), and the strength of 
the interaction, usually parameterized by {kpa)~^. Throughout this work we shall use units in which Boltzmann's 

OO . constant is fc^ = 1. The function ^ represents the ratio of the energy E to the energy of a free Fermi gas at the same 

■ density Ep — 3/5iNep. The value of ^ at unitarity and at T = 0, which we shall denote ^s, has been determined 
^ , approximately by various authors, and recent Quantum Monte Carlo calculations and extrapolations to zero range 

r*"! yield = 0.40(1) (see and references therein). Besides ultracold atomic gases, the unitar y re gime is relevant for 
rN ] dilute neutron matter, although in that case finite effective range effects cannot be neglected |15| . 

■ In this paper the determination of the universal function £^{T/ep, l/kpo), along with other thermodynamic quanti- 
■ ■ ' ties (including the critical temperature Tc for the onset of superfluidity), will constitute our main results. We explore 

the unitary limit, where kp\a\ —> oo, as well as the case of finite scattering length in the range —0.5 < {kpa)~^ < 0.2. 
The paper is organized as follows: in section II we formulate the problem and describe the nonperturbative technique 
based on the discretization of the space-time and subsequent evaluation of the thermodynamic quantities through 
a Quantum Monte Carlo simulation. In section III we describe the numerical and computational techniques that 
were used. The main results are discussed in sections IV (at unitarity) and V (away from unitarity), and the final 
conclusions in section VI. 
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II. MATHEMATICAL FORMULATION: FERMIONS ON A SPACE-TIME LATTICE 



A. The Hamiltonian and the running coupling constant 



The interaction that captures the physics of a dilute, unpolarized system of fermions is a zero-range two-body 
interaction V(ri — — —gS(ri — r2). The Hamihonian of the system in the second quantization representation 
reads: 



H 



^2y2 



(ii.i) 



where ^^(r) = ^j;'^ {r)'tjj^{r). This interaction requires that we specify a regularization procedure, which we do by 
introducing a momentum cut-off hkc (thus requiring all two-body matrix elements to vanish if the relative momentum 
of the incoming particles exceeds the cut-off). Once this cut-off is imposed, the value of the bare coupling g can be 
tuned to fix the value of the physical renormalized coupling, which in this case will be the s-wave scattering length a. 
Indeed, the diagonal T matrix describing two-particle scattering induced by the interaction takes the simple form: 



T(fc) = 



'27r)3 ^ 









(27r)3 \ 



(27r)3 



1 + 



gm 



kc 
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(n.2) 

which is equivalent to finding the vacuum 4-point amplitude and determining the scattering length by summing all 
the 'bubble' diagrams (see Fig. [1]). The low-momentum expansion of the scattering amplitude reads: 



—ik 



gm 



2k c 



2fc2 



n -1 



Oik') 



(II.3) 



At low momentum we have, by definition, f{k) = [—ik — l/a + rcsk'^ /2 + 0{k^)] ^, which gives the relation between 
the bare coupling constant g and the scattering length a at a given momentum cutoff hkc- 



m 



+ 



kcvn 



(11.4) 



Note that an effective range is generated which is independent of the coupling constant r^f — — — . 

nkr 



XO-OC 



FIG. 1: n-th term in the bubble sum. 



B. Discrete variable representation 

To determine the thermal properties of spin 1/2 fermions in a non-perturbative manner, we have placed the system 
on a three dimensional (3D) cubic spatial lattice with periodic boundary conditions. The lattice spacing I and 
size L = Ngl introduce natural ultraviolet (UV) and infrared (IR) momentum cut-offs given by hkc — T^h/l and 
HAq = 27rfi./L, respectively. The momentum space has the shape of a cubic lattice, with size 2-Kh/l and spacing 
2nh/L. To simplify the analysis, however, we place a spherically symmetric UV cut-off, including only momenta 
satisfying k < kc < n/l. 
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The discretization transforms the functions of continuous variables into functions defined on a discrete set of 
coordinate values: 

Via(r.) -> ^^(i) (II.5) 

i^+(rO ^ V'+li) (11.6) 

flaiVi) no-(i), (n.7) 

where = il and i = (ixjij/j^z) denotes the lattice sites, and ixiiy,iz — l,---,-^s- The discretization can affect 
the accuracy of the obtained results and therefore requires careful analysis. We address this issue by discussing the 
so-called discrete variable representation (DVR) basis sets, which is the underlying framework of our lattice approach, 
see Ref. [H]. 

Let us call Ti = L?'{M) the Hilbert space of our problem, i.e. the set of square-integrable wave functions on the 
manifold M . Let P be a projector defined on Ti such that S = PTi, is the projected subspace. Given a set of N grid 
points {xq, a = 0, N— 1} in d dimensions one can define projected (5-functions: Aa{x) = P[S{x~Xa)]- Alternatively 
|Aq) = P\xa) or Aa{x) = {x\Aa), in Dirac notation. 

It follows that 

{A„\Af3) = Afj{x^) = Alixfj) (IL8) 
and as a result the set of projected (5-functions {Aq,, a — 0, — 1} is orthogonal if and only if 

Ac,{xp) = Kc,6c.(3, {11.9) 
where — (Aq|Ac), and we can normalize the projected functions to get their orthonormalized version: 

|F,) = ^^|A„). (ILIO) 



Given a wave-function -0 g 5, an expansion of the form 



V'(x)= Vc„F,(x) (ILll) 



exists, and the coefficients are given by the values of ip{x) at the lattice sites Xa- 

dxF*{x)^{x) = -l=^{x^) (IL12) 

V Aa 

On the other hand, if il){x) is not fully contained in S, then our basis set will not be sufficiently rich to represent 
it. However, if the semiclassical region of phase space that we wish to represent (see Fig. [5]) is contained in S (and in 
particular if the UV momentum cutoff hkc is larger than the highest momentum we wish to represent, or equivalently 
if the lattice spacing / is chosen to be sufhciently small), then it can be shown (see Ref. [ia|) that the errors are 
exponentially suppressed. 

The representation of the Hamiltonian on the lattice has been obtained by noticing that two terms representing 
the kinetic and interaction energy are local in momentum and coordinate spaces, respectively. Thus, the kinetic term 
reads: 

^-^ 2m 

'y=\A T 

where i = {i^, iy, iz) enumerates the lattice sites in the momentum space, ix, iy, iz = 0, ±1, ±2, Ns/2(ot — Ns/2) and 
kr — VItt/ L. The operator na{T^) denotes the occupation number operator of the single particle state with momentum 
?ikj and spin <t. On the other hand, the interaction becomes a simple Hubbard attractive potential 

F = -g/3^nT(iK(i)- (H.M) 

i 

From this point on we shall omit any factors of or coming from the volume elements in real and momentum 
space, which amounts to a specific choice of units. If the cutoff is chosen to be fee = tt/^, then infinite scattering length 
corresponds to a couphng given hy g = 2'Kh^l/m (see Eq. (|IL4p ). 

In order to avoid numerical inaccuracies associated with the discretization of the differential operators on the lattice, 
we use both momentum and coordinate representation of the lattice and a Fast Fourier Transform (FFT) to switch 
between the two. In the following we shall rename the indices i and i in favor of the more suggestive names k and r, 
respectively. 
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FIG. 2: (Color online) Representation of phase space. The dashed line shows a classical trajectory, while the shaded region 
represents the quantum fluctuations around such trajectory. The horizontal and vertical lines show the cutoffs in real space L, 
and in momentum space hkc = hn/l, where I is the lattice spacing (assuming the lattice is a uniform square lattice). Functions 
within the larger oval region are well described by a given basis if the cutoffs enclose that region, as shown in this figure. 



C. Discrete auxiliary fields and positivity of the probability measure 

In order to study the thermal properties we chose the grand canonical ensemble, where the thermodynamic variables 
are the temperature T, the chemical potential fi and the volume V . The partition function and the average of an 
observable O are calculated according to 



Z(/3,^,y) - Tr{exp[-/3(i/-AiiV)]}, 
Tr |dexp[-/3(i/-/i7V)]| 



0{p,^l,V) 



(11.15) 



where /3 = 1/T (as mentioned before, in this work we will take Boltzmann's constant to be fc^ = 1). By factorizing 
the statistical weight using the Trotter formula, one obtains 



exp[-/3(iJ - fiN)] = Yl cxp[-T{H - fiN)] 



(11.16) 



where (3 = NtT. The next step is to decompose the exponentials on the right hand side into exponentials that 
depend on the kinetic and potential energy operators separately. This can be achieved to a suitable order through 
the factorization 



M 



exp[-T(i + B)] ~ exp[-afcTi] exp[-6fcTi?], 



(11.17) 



where the coefficients Uk , bk are determined by the required order of accuracy fv?\ . However in order to obtain stable 
numerical results all of the coefficients have to be positive. This poses a practical limitation on the above formula 
which works only up to the second order, for which the expansion is 



exp[~T(iJ — fiN)] = exp 



t{K - fiN) 



exp(— tV) exp 



(11.18) 



where K is the kinetic energy operator, whose dispersion relation, for momenta smaller than the cut-off, is given by 
£k — h'^k'^/2m. It is important to note that, because we have used the expansion up to C(t^), when calculating the 
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partition function this becomes ©(t^). Indeed, the statistical weight involves a product of A'^,- factors and is given by 
the following expression: 



exp[-/3(i7 - fiN)] = exp 



t{K - ^^N) 



exp[— tF] exp[— t(X — ^7V)] | exp 



-{K - fiN) 



+ 0{t^) (11.19) 



Notice also that this approach docs not depend on the choice of dispersion relation. One may choose to implement a 
discrete derivative for the kinetic energy (which corresponds to the Hubbard model) based on the second difference 
formula: 



f{x + l) + f{x-l)-2f{x) 
P 



(11.20) 



This results in a dispersion law e^. ~ siv?{kxl/2) + siv?{kyl/2) + sin^(fc2//2), where I is the lattice spacing. Notice 
that this dispersion relation is not spherically symmetric in momentum space. In that case the energy of the higher 
momentum states differs noticeably from the continuum case. Indeed, as we show in Fig. [3l the dispersion relation 




FIG. 3: (Color online) The solid blue line shows the dispersion relation used in this work and the dashed lines and purple area 
result from a lowest order second difference discrete derivative (see text for discussion). The units in this plot are set by the 
lattice spacing: the wavevector k is in units oi l/l and the energy in units of h^/2ml^. 



that results from the second difference formula deviates significantly from behavior already at fc ~ 7r/2Z and as a 
result the number of physically meaningful available states is only ~ 1/2^ of the whole phase space. Furthermore, the 
deviation from the fc^ law depends on the angle in momentum space, and as a result the discrete derivative formula 
sweeps the shaded area in Fig. [3l This is not the case with our choice of ek as we shall consider the kinetic energy 
operator in momentum space. According to our experience this indeed minimizes the discretization (high momentum) 
errors. 

The interaction factor can be represented using a discrete Hubbard-Stratonovich transformation, similar to the one 
in Ref. [H: 



exp[gTn|(r)n|(r)] 



E 

cr(r,Tj)=±l 



[1 + Aa{r,Tj)h^{r)][l + A<j{r,Tj)hi{r)]. 



(11.21) 



where A = y/ exp((7r) — 1, tj labels the location on the imaginary time axis, where j = 1, ...,7Vr, and a{r,Tj) is a 
field that can take values ±1 at each point on the spacetime lattice (the name of this field should not to be confused 
with the spin variable, conventionally also called a). This identity can be proven simply by evaluating both sides at 
n{l l}(r) = 0, 1. This discrete Hubbard-Stratonovich transformation is sensible only for A < 1, which means that the 
imaginary time step which cannot exceed log 2. In practice however the actual value of r has to be much smaller 
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due to the finite size of the r step in formula pi.l9[) (see also next section) . The fact that a takes only discrete values 
is a consequence of Fermi-Dirac statistics, and results in a discrete configuration space for the field a. This is the 
main advantage for numerical applications as compared to the case of a continuous Hubbard-Stratonovich field [l8| . 
Taking all this into account, the grand canonical partition function becomes 



Z(/3,^,y) = Tr{exp[-/3(i7-M^)]}= j \{Vcj{v,Tj)TY lJ[{a}). 



(11.22) 



where we define 



and 



W,(M) = exp 



ma})^\{W,{{a}) 



+ Aair, Tj)h^ (r)] [1 + Aa{T, Tj)h^ (r)] exp 



t{K - nN) 



Since a is discrete, the integration is in fact a summation: 



{T(r,Ti)}=±l{CT(r,T2)}=±l {(T(r,rjvJ} = ±l 



where 



E = E E - E ' 

{(T(r,rj)} = ±l ;7((l,0,0),r3) = ±l a({2,0,0),rj)=±l a{{N^ ,Ns ,N^) ,Tj)=±l 

In a shorthand notation we will write 

U{{a}) ^ T^expj- / dT[h{{cr}) ~ fiN] 



(11.23) 



(11.24) 



(11.25) 



(11.26) 



where T,- stands for an imaginary time ordering operator and h{{a}) is a resulting a-dependent one-body Hamiltonian. 
It is crucial to note that U{{a}) can be expressed as a product of two operators which describe the imaginary time 
evolution of the spin-up and spin-down fermions: 



(11.27) 



The operators for spin-up and spin-down are identical for the case in which the chemical potential is the same for 
each spin: /Uj = = /x (which is the case we consider in this work) . 
The expectation values of operators take the form 



0(/3,M,^") 



Tr [dexp[-P{H^fiN)]} ^ /• I];, Pa(r, )Tr Zi(M) Tt du{{a}) 



Z(/3,/i,V) 



TTU{{a}) 



(11.28) 



where we have introduced Tr U{{a}) for convenience: in the numerator it represents the probability measure used in 
our simulations (see below), and in the denominator it serves the purpose of moderating the variations of Tr OU{{(7}) 
as a function of the auxiliary field a. 

All of the above traces over Fock space acquire very simple forms [l^, [23] (see next section), and can be easily 
evaluated. In particular, TrU{{a}) can be written as 



TrU{{(T}) = det[l -t-Z^({CT})] = det[l +Ui{{a})] det[l -t-Z^|({cr})], 



(11.29) 



where 14 (without the hat) is the representation of U in the single-particle Hilbert space. The second equality is a 
result of the decomposition pi.27p . This identity is easy to prove by expanding both sides, taking into account that 
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W is a product of exponentials of 1-body operators. In the case considered in this work the chemical potential is 
the same for spin-up and spin-down fermions, so it follows that det[l = det[l -|- Z//|({cr})]. This implies 

that Tr W({(t}) is positive, i.e., that there is no fermion sign problem for this system. Indeed, this allows to define a 
positive definite probability measure: 



TrUiM) {det[l+^^T(W)]}^ 
Z{f3,fi,V) 



1 



Zif3,^i,V) 



exp(2 tr (log[l+^(M)])) 



(11.30) 



where the exponent in the last equation defines the negative of the so-called effective action. 

The many-fermion problem is thus reduced to an auxiliary field Quantum Monte Carlo problem, to which the 
standard Metropolis algorithm can be applied, using Eq. (jIL30|) as a probability measure. Before moving on to 
the details of our Monte Carlo algorithm, we briefly discuss the expressions used to compute a few specific thermal 
averages. 



D. Calculation of observables 



Let us consider the one body operator 



0= J2 / d'rid3r2V'+(ri)at(ri,r2)Vt(r2) 



(11.31) 



From Eq. (|TL28)) it follows that 



(6).i:p,w,^fM).i:p,w,.^'»<w) 



TrUiia}) 



' det[l+ U{{a})Y 



The calculation of the last term requires the evaluation of 



Tr 



^+(ri)^t(r2)Zi(W) = Sstdet[l +U{{<j})]^ns{r^,r2,M) 



where s,t denote spin (f or j), and 



,(ri,r2,{cr}) = Vki(ri) 



1+Usi{a}) 



(11.32) 



(11.33) 



(11.34) 



ki,k2 



Here (/3k(r) = exp(ik • r)/L^/^ are the single-particle orbitals on the lattice with periodic boundary conditions. This 
holds for any 1-body operator O, \i U is a product of exponentials of 1-body operators, as is the case once the 
Hubbard-Stratonovich transformation is performed. It is then obvious that the momentum representation of the 
one-body density matrix has the form 



,(ki,k2,{cr}) = 



Us{{a}) 



l+Us{{cj}) 



(11.35) 



ki ,k2 



which, for a noninteracting homogeneous Fermi gas, is diagonal and equal to the occupation number probability 

l/(exp[/?(ek — m)] + 1) of a state with the energy = — • 

Summarizing, the expectation value of any one-body operator may be calculated by summing over samples of the 
auxiliary field (T(r,Tj): 



(O) = j \{Vcj{v,T,)P{{a}) J2 E a.(ri,r2)n,(ri,r2,{a}) 



(11.36) 



ri,r2 8=1,1 

In particular the kinetic energy can be calculated according to 
rirr, ^''^(i-, ^.)Tr Uila}) Tr XZ^({a}) 



(K) 



k<ka 



l[Vair,T,)P{{a})Y^ ^ 

k s=T,i 



ns(k, k, {a-}) 



2m 



(11.37) 
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Analogously, for a generic two-body operator: 

J c?'r;d3r^d3rid3r2V'+(ri)V't+(r^)0«t™K,r^,ri,r2)^„(r2)^,(ri) 

s,t,u,v=l,'\ 

in order to calculate (O) one needs to evaluate expression 
Tr [^+(r'i)V3+(r'2)Vi„(r2)^„(ri)Zi({a}) = 

= (det[l +U{{<T})]f {5suStvns{r[,ri, {cr})nt{r'2,r2, {a}) - SsvStuns{r[,r2, {(j})nt{r'2,ri, {a})) 
Hence, for the expectation value of the two body operator we get 



{d)= fl[Va{r,T,)Pi{a}) ^ ^ 



Ostst(ri, ri, r2)n5(ri, ri, {cr})nt(r2, ra, {ct}) 
Ostts (r'l, r'a, ri, r2)n^(r'i, ra, {cr})nt(r'2, ri, {ct}) 



In particular, the expectation value of the interaction energy reads: 

{V)=-g j n2?'^(r,r,)F(W)^nt(r,r,{a}K(r,r,M) 

It should be noted that in the symmetric system = /ij) 

n]{Yy,{a}) = n|(r,r', {cr}). 

Hence, 



{V) - -3 / n ^^(^' ^.O^Ci'^}) E["T(r, r, {a})] 



It is useful to introduce the correlation function 

2 



92{v) = (^^) y'd3rid3r2(^t(j.^+r)Vl(r2+r)^i(r2)^T(ri)) 
X j (f'Yid^V2n^{vx +r, ri,{cr})nx(r2 +r,r2,{cr}) 



(H.38) 



(H.39) 



(H.40) 



(11.41) 



(11.42) 



(H.43) 



(11.44) 



(where N is the average particle number) which is normalized in such a way that for a noninteracting homogeneous 
Fermi gas 52 (r) — S '^^, ^ — and 52(0) — 1. 



III. NUMERICAL METHODS AND COMPUTATIONAL ISSUES 



A. Metropolis Monte Carlo 



Once we have written the observables as in Eq. (jll.28p . the next step is to sum over all possible configurations of 
cr(r, Tj). For a lattice size x Nr (where typically = 8 and Nr ~ 1000), performing the sum over the 2^^^^^ 
points in configuration space, is an impossible task for all practical purposes. It is in these cases that a Monte Carlo 
(MC) approach becomes essential. By generating M independent samples of the field a{r,Tj) with probability given 
by pi.30[) . and adding up the values of the integrand at those samples, one can estimate averages of observables with 
0{1/VJV) accuracy. 
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The standard Metropolis algorithm was chosen to generate the samples. At every MC step the sign of a was changed 
at randomly chosen locations on the space-time lattice. An adaptive routine increased or decreased the fraction of 
sites where a was updated, so as to maintain an average acceptance rate (over ~ 100 consecutive Metropolis steps) 
between 0.4 and 0.6. 

In order to compute the probability of a given a configuration, it is necessary to find the matrix elements of 
which entails applying it to a complete set of single-particle wave- functions. For the latter we chose plane waves (with 
momenta Kk < hkc). This choice is particularly convenient because one can compute the overlap of any given function 
with the whole basis of plane waves by performing a single Fast Fourier Transform (FFT) on that function. 

In practice, the calculation of the matrix elements of U proceeds by evolving the above mentioned set of wave- 
functions in imaginary time, applying the operators Wj in Eq. pi.23p in sequence. The operator Wj is in turn 
made up of kinetic energy and potential energy factors (see Eq. (|II.24p ). in a specific order. The application of such 
factors was implemented in momentum space and real space, respectively, using FFT to switch between them. 



B. Singular value decomposition 



80r 
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FIG. 4: Left Panel: Condition number of W as a function of f3. Squares: with SVD. Crosses: without SVD. Right Panel: 
Convergence of simulated particle number relative to exact solution A'^o, as a function of time step r in units of to, defined 
in the text. 



Matrix multiplication, especially in the form of FFT's, is ubiquitous in our algorithm. It is well known that matrix 
multiplication is numerically unstable when the matrices involved have elements that vary over a large range. This 
is true in our case at low T, with exponentially diverging scales. To avoid instabilities it is necessary to separate the 
scales when multiplying matrices, and the Singular Value Decomposition (SVD) technique serves such purpose. In 
this section we follow the same approach developed in Ref. 19] to introduce SVD in our calculations. 

Let us write the matrix U{{a}) more explicitly: 

^iW}) = n = ^N^^N^^i-mwi (111.1) 

where the Wfc({cr}) will be x matrices, for a single-particle basis of dimension N. Let us then define 

Ua = 1 (III.2) 
Ui = Wi 
U2 = W2W1 
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To separate the scales one decomposes the matrix Un-i before muhiplying it by Wn to get Un- This process begins 
as foUows 

Uo = 1 (III.3) 

U2 = W2W1 = (W25iPi)Vi = S2V2V2V1 

where Si and Vi are orthogonal matrices (not necessarily inverse of each other), and I?i is a diagonal positive matrix 
containing the singular values of Ui. The idea is that the actual multiplication should be done by first computing 
the factor in parenthesis in the last equation. This factor is then decomposed into S2'D2V2i m preparation for the 
multiplication by Wa, and so on. A generic step in this process looks like this: 

U,, = WnUn-l = W„5„_iI?„_iV„_iV„_2...Vl (III.4) 

so in the end 



= U{{a}) = SnVn^Vn^Vn^-i-Vi = SVV (III.5) 

where we have decomposed the full product in the last step. Calculation of the determinant, and therefore of the 
probability measure, becomes straightforward if we perform one more SVD, in the following chain of identities: 

det(l +U{{cf})) = det(l + SVV) = det(5(5t + V)V) ^ det(5 SVV V) = det(55) det(I?) det(VV) (III.6) 

For equal densities (symmetric case), where we need this determinant squared, we only care about the factor in the 
middle of the last expression, since the other two are equal to 1 in magnitude. Indeed, in that case we can write the 
probability measure as 

P(M) = exp|^^logJ,j (IIL7) 

where di > are the elements in the diagonal of V, and M is the dimension of the single particle Hilbert space. 
The SVD decomposition is useful when evaluating occupation numbers. Indeed, 

"^^"^^^ = 1 = 1 - vtytp-i^t^t (ni.s) 



l+U{{a}) l+U{{a}) 

which is very easy to compute since the matrix V is diagonal. 

In Fig. [4] (left panel) we show the behavior of the condition number (defined as the ratio of largest to smallest 
eigenvalue) of the matrix Un^ in the single-particle Hilbert space. The number of SVD's required to stabilize the 
calculation grows as we increase (3. If no SVD's are used, the condition number saturates, indicating loss of information 
due to poor separation of scales in matrix multiplication. In our calculations we have made limited use of SVD, ranging 
from 2 decompositions at the highest T to 8 decompositions at low T's. 



C. Tests and cross-checks 



In order to verify the correctness of our code we performed several tests. As a first check the thermodynamics of a 
free gas was reproduced when setting g = 0. This is an elementary test, as in this case the MC part of the algorithm 
is superfluous. 

To check our results at g 7^ we diagonalized the Hamiltonian exactly, restricting the phase space to the lowest 
7 single-particle momentum states. This entails constructing all of the 2^^ states [2^ for each spin) and computing 
and block diagonalizing the Hamiltonian (which comes in blocks identified by fixed particle numbers (iV|,iV|)). An 
average desktop computer can complete the whole task in about 5 minutes. This test provided an estimate for the 
size of the step in the imaginary time direction. In Fig. [4] (right panel) we plot the difference between the simulated 
number of particles N and the exact value iVo, as a function of the time step, in units of tq = In 2/(7, with all the 
other parameters (T, /i, etc) fixed. We conclude from this data that it is safe to take r = tq/IO, as the relative error 
falls below 10%. 
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In the unitary limit universality implies thermodynamic relations that provide self-consistency checks. For instance, 
it can be proven that ii E = 3/5epN(_{T/ep) then 

^^d^)~'-^r^dy^-^. (III.9) 

In our simulations /i, T and V are input parameters, while E and sp are part of the output. We have checked that 
the above relation is satisfied by our data. 



D. Density dependence and the role of periodic boundary conditions and the high- momentum cutoff 



An MC simulation of a Fermi gas in the unitary regime makes sense only if the following conditions are satisfied: 

27r 

T 



Ao 



<^kp -^kc 



n 
1' 



(III.IO) 



The size of the box L defines the lowest momentum scale HAq — 2'Kh/ L one can resolve in such a simulation, while 
the lattice constant I defines the smallest inter-particle separation accessible. To better appreciate how the restriction 
(|IIL10[1 affects the simulations we present in Fig. [5] the errors one incurs due to both the upper momentum cutoff, 
needed to regularize the two-body interactions, and the use of boundary conditions. We calculate the total particle 
number in a box with periodic boundary conditions in three different ways: 



k<ka 
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1 - 


e{k) 


+ U-fi 
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h A2 



(III.ll) 

(III.12) 
(III. 13) 



where e(fc) is the single-particle kinetic energy, k = 2Tr /L{nx,ny,nz) and rix^y^z are integers, and using parameters 
characteristic for the unitary Fermi gas at zero temperature, see Eq. (jV.4[) and Ref. [2l| . Thus N is the exact particle 
number in such a box if no restrictions on momenta were imposed, Ncont would be the actual particle number only 
if all momenta smaller than the upper cutoff momentum hkc = ttH/I would be taken into account, and Niatt is the 
particle number one would obtain if the simulation were performed in a box with periodic boundary conditions and an 
upper cutoff momentum hkc- The quantity | — Ncont\/N is thus a measure of the error introduced by a finite cutoff 
momentum hkc alone. It is debatable whether one has to take into account in Ncont the contributions from momenta 
larger than the cutoff momentum hkc, since the physics at those momenta is not simulated correctly anyway. In a 
box with periodic boundary conditions by default one can access only a finite set of discrete momenta. As it is clear 
from the figure, imposing periodic boundary conditions alone leads to very large errors especially in the case of small 
densities, when the Fermi momentum hkp = h{i'K'^N / L^y-^^ becomes smaller than the momentum ?iAo = 2nh/L. It 
is important to notice as well that the magnitude of these errors are independent of the presence or absence of the 
upper cutoff momentum hkc- A rule of thumb would suggest that one has to have at least ten particles in a box 
in order to keep the errors at an acceptable level. This is relatively easy to understand physically, as in a box with 
periodic boundary conditions the lowest single-particle levels are: one level with zero energy, six levels with energy 
h?K'^/2m., and twelve levels with energy 2h^K^/2m, not counting spin degeneracy. 

The quantity kpL/2TT is basically the ratio between half of the box size (£/2) and the average inter-particle 
separation (w n/kp), and in our simulations we were able to probe kpL/2TT « 2.5. In Ref. 5^ this ratio was increased 
to about 5, i.e. the system was significantly more dilute. Notice however, the errors incurred in a calculation in a box 
with periodic boundary conditions and very small particle number (and thus very low density) become unacceptably 
large, as demonstrated in Fig. (5] 

With the densities used in this work, we find excellent agreement between our results and GFMC and DMC results 
0) H) 0) IE [Hi [3 for essentially all of the quantities of interest. All these calculations however were performed 
typically for particle numbers between 10 0, y| and less than 66 0, H, |3, H, [H, [i3l , where the errors arising from 
imposing periodic boundary conditions are minimal. Only the long range universal critical behavior of the two-body 
density matrix discussed in Section IV {G{r) cx r~'^+''') requires very low filling factors, but in very large box sizes, 
and in a very narrow temperature interval around the transition. 
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FIG. 5: . (Color online) The error in particle number obtained in a simulation in a box with periodic boundary conditions and 
upper momentum cutoff \Niatt — N\/N, solid line (blue), compared with the error obtained in a truly continuum model with 
no periodic boundary conditions imposed but only an upper momentum cutoff \Ncont — N\/N, dash-dotted line (red), and the 
"numerical" apparent error \Niatt — Ncont\/N, dotted line (green). We display here two cases, of simulations in a box with 
lattice sizes 10^ and with 20^. In the case of the quantity \Niatt — Ncont\/N we see no lattice size dependence. 



IV. RESULTS IN THE UNITARY LIMIT 



The results of our simulations for lattices ranging from 



1732 (at low T) to 8^ x 257 (at high T) and for 



2 • • • 20 X 10^ Monte Carlo samples (after thermalization) are shown in Fig.[71 The imaginary time step was chosen as 
T = min(m/^/157r^?i^, In 2/10(7), The first bound comes from the inverse of the highest single-particle kinetic energy 
available on the lattice, namely Ek^max = 3/i^7r^/2m/^, and the second bound results from our specific choice for 
the discrete Hubbard-Stratonovich transformation (in both cases an extra factor of 1/10 was included based on our 
comparison with an exact calculation, as described in Section fill CI ) The Monte Carlo autocorrelation length was 
estimated (by computing the autocorrelation function of the total energy) to be approximately 200 Metropolis steps 
at T « 0.2eF- Therefore, the statistical errors are of the order of the size of the symbols in the figure. The chemical 
potential was chosen so as to have a total of about 45 particles for the 8^ lattice. We have performed however 
calculations for particle numbers ranging from 30 to 80, for lattice sizes 8'^ and 10^ and various temperatures, and in 
all cases the results agreed within discussed above errors and systematics. In all runs the single-particle occupation 
probabilities of the highest energy states were well below a percent for all temperatures. This can be seen in Fig. [SI 
where the dispersion in the data at fixed momentum is the result of statistical errors and the fact that the lattice is 
not spherically symmetric. 

We tried to extract from our data the asymptotic behavior in the limit of large momenta n{k) cx C{kp /k)" , which 
according to the theory [H, should at all temperatures be governed by the same exponent, namely n = 4. Our 
results are consistent with a value of the exponent n = 4.5(5), thus in reasonable agreement with the theoretical 
expectation. 

All the quantities computed present a number of common features that are easily identified, in particular a low 
and a high temperature regime, separated by a characteristic temperature that we estimate to be Tg = 0.23(2)£i?. 
We shall discuss in the next sections whether Tq can be interpreted as the critical temperature for the onset of 
superfiuidity. 
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FIG. 6: (Color online) Occupation probability in y-log scale at intermediate temperatures, T — 0.20eF, in blue circles for 8'^ 
and red squares for 10'^. 



A. Energy and chemical potential 

As T ^ the energy tends to the T = results obtained by other groups [2, [3, Q • This confirms such results, as 
the algorithms they used (namely Green Function Monte Carlo) are constrained by the existence of a sign problem, 
which is not the case in the present approach. The solid line in Fig. [7] shows the energy of a free Fermi gas, with a 
shift down given by 1 — with ^„ = + (5^ ~ 0.55, where S,s = 0.4 and 



5E 



A 



0.15 



(IV.l) 



is the condensation energy in units of the free gas ground state energy. Here we have used the BCS result 5E = 
(see Ref. [HI), where A is the pairing gap found in Ref. [l^ to be A ~ 0.50eF- One can also find the value of ^„ from 
our data by determining what shift is necessary to make the solid curve in Fig. [7] (which corresponds to the free gas) 
coincide with the high temperature data (where the gas is expected to become normal). We find that such procedure 
gives ^„ ~ 0.52, which is roughly consistent with the value quoted above. This number should also be compared with 
the results of Refs. 0,111], namely ^ 0.54, and Ref. [131, that finds ^„ ~ 0.56. 

For T < To we observe that the temperature dependence of the energy can be accounted for by the elementary 
excitations present in the system in the superfluid phase: boson-like Bogoliubov-Anderson phonons and fermion-like 
gapped Bogoliubov quasiparticles. Their contribution is given by 



Eph+qp{T) — '^^pN 




(IV.2) 



(IV.3) 



where A is the approximate value of pairing gap at T = determined in Ref. Q to be very close to the weak-coupling 
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FIG. 7: Left panel: total energy E{T) with open circles, and the chemical potential jJ.{T) with squares, both for the case of 
an 8^ lattice. The combined Bogoliubov- Anderson phonon and fermion quasiparticle contributions Eph+qp{T) (Eq. (IIV.2|l 'l is 
shown as a dashed line. The solid line represents the energy of a free Fermi gas, with an offset (see text). Right panel: entropy 
per particle with circles for 8^ lattice, and with a dashed line the entropy of a free Fermi gas with a slight vertical offset. The 
statistical errors are the size of the symbol or smaller. 



prediction of Gorkov and Melik-Barkhudarov [2^, and ~ 0.44 Q and ep — h?kp/2m and n = kp/Sir^ respectively. 
Notice that the estimate for central value of has decreased in the last years: Ref. [l^ reported = 0.42(2), 
while recently = 0.40(1) has been reported in Ref. [3, H^, and = 0.37(5) in this work (see Table U below), 
even though all these results agree within quoted errors. The sum of the contributions from the excitations, namely 
Eq. IIV.21 is plotted in Fig. [7] as a dashed line. Both of these contributions are comparable in magnitude over most 
of the temperature interval (To/2,ro). Since the above expressions are only approximate formulas for T <C Tc, the 
agreement with our numerical results may be coincidental. 

At T > Tc the system is expected to become normal. If Tq and Tc are identified, then the fact that the specific heat 
is essentially that of a normal Fermi liquid Ep(T) above Tq is somewhat of a surprise, as one would expect the presence 
of a large fraction of non-condcnscd unbroken pairs. Indeed, the pair-breaking temperature has been estimated to be 
T* ~ 0.55£f, based on fluctuations around the mean-field, see Refs. [sOjUH, implying that for Tc < T < T* there 
should be a noticeable fraction of non-condensed pairs. As we shall see, Tc ~ O.l^ep < Tq (this result for Tc was first 
obtained in Ref.Q), and as one can see in the caloric curve of Fig. [3 the specific heat deviates from the normal Fermi 
gas, in the range Tc < T < Tq. 

On the other hand, the chemical potential fi is almost constant for T < Tq, a fact reminiscent of the behavior 
of an ideal Bose gas in the condensed phase, even though in such a phase our system is strongly interacting and 
superfluid. This implies a strong suppression of fermionic degrees of freedom at those temperatures. Moreover, 
assuming /i(r) = const, for T < Tq implies that 

E(r, = 4.,f(i:), f(Z-).f,+c(^)". (IV.4) 

which is the temperature dependence of an ideal Bose condensed gas. According to our QMC results, the value of n 
extracted from our data is n — 2.50(25). 



B. Entropy 

From the data for the energy E and chemical potential fi one can compute the entropy S, because in the unitary 
limit the relation PV = holds (true of a free gas as well), where P is the pressure, V is the volume and E is the 
energy. It is straightforward to show that 



S _E + PV - iiN _ ^(x) - C(a;) 
iV ~ WT ^ X 



(IV.5) 
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FIG. 8: Entropy per particle as a function of the energy per particle ^. The dotted line shows the location of Tc (see next two 
sections); the dashed line shows the location of Tq. 



where = fx/ep and x = T jsp which determines the entropy per particle in terms of quantities we know from our 
simulation. The entropy also departs from the free gas behavior as the temperature is lowered below Tq. 



As indicated in 
suggestion of Ref. 



, this data can be used to calibrate the temperature scale at unitarity. Indeed, extending the 
I, from known T in the BCS limit, the corresponding S{Tbcs) can be determined. Then, by 
adiabatically tuning the system to the unitary regime, one uses S{Tbcs) — SlTunitary) to determine T at unitarity (In 
practice the experimental procedure goes in the opposite direction, namely measurements are performed at unitarity 
and then the system is tuned to the deep BCS side, see Ref. [11]). 

On the other hand, knowledge of the chemical potential as a function of temperature (see previous section) allows 
for the construction of density profiles by using of the local density approximation. In turn, this makes it possible 
to determine S{E) for the system in a trap, fixing the temperature scale via dS/dE — 1/T. Direct comparison with 
experimental results, in remarkable agreement with our data, has been demonstrated by us in Ref. 35]. In Fig. [8] 
we show the entropy per particle S/N as a function of the energy per particle ^ for the homogeneous system (See 
Ref. [11] for the corresponding result for the case of the unitary gas in a trap.) 



C. Two-body density matrix and condensate fraction 



Information about the onset of critical behavior (e.g. a superfluid phase transition) can be obtained by studying 
an appropriate order parameter, both as a function of temperature and system size. In the case of superfluidity in 
two-component Fermi systems, which is a particular example of off-diagonal long-range order [36j . the order parameter 
is the long-distance behavior of the two-body density matrix g2{r), defined in Eq. pi.44[) . At unitarity, knowledge 
of g2{f) is enough to determine the condensate fraction a = lim^^oo y52('')- On the BCS side of the resonance, as 
discussed in [^ , the calculation of a demands also knowledge of the one-body density matrix: 

p{r) = ^ y d'ri(V|(ri + r)VT(ri)) (IV.6) 

and a = lim^^co t52('') - 
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FIG. 9: (Color online) Left Panel: Projected two-body density matrix (see text) as a function of position on the lattice. In 
solid blue for a 10^ lattice, for the temperature range (O.Ief, O.Ssf). The T = results from Ref. 5] are shown with green 
circles and error bars. Right Panel: Condensate fraction a(T), black for 6'' (highest), red for 8^ (middle) and blue for 10"^ 
(lowest) lattices, respectively. 



In Fig. [9] we show our results for g2{r), as a function of the dimensionless lattice position kpr (left panel) and the 
extracted condensate fraction a at unitarity for several lattice sizes (right panel). In Fig. I10[ on the other hand, 
we show schematically the generic form of a correlation function G(r/^corr)(such as g2{r)). At short distances, r ~ I 
(where here / should be regarded as an intrinsic short distance scale of the problem) , the behavior of G depends on 
the system under consideration, i.e., it is non- universal. In the region I <^ r < ^corr (which extends to infinity at a 
phase transition because there ^corr ~^ oo, see next section) the form of G is universal, in the sense that it depends, 
quite generally, only on the spatial dimensionality of the problem and the internal symmetries of the Hamiltonian. In 
that region G - r-^'^+v)^ where r] ~ 0.038 is a universal critical critical exponent (see e.g. [37j ). Finally, for r > ^corr 
the correlation decays exponentially with r/^corr- 
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FIG. 10: (Color online) Generic form of a correlation function G(r/^c, 
the correlation length ^corr. 



as a function of the radial coordinate r in units of 



Below the critical temperature, where a is non- zero, its value can be extracted from the asymptotic form of g2{f)- 
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In Fig. [9] (left panel), one can see that the asymptotic value of (72(?') for 10^ is approached even for r < L/2, so 
the lattice sizes used are reasonable large to determine a. It is possible that for larger lattice sizes a could become 
smaller. On the other hand, there is virtually no room for finding the power law that characterizes the universal 
critical behavior of this function (see Fig. [TO]). In other words, corrections to universal scaling will be important close 
to the critical point for these small volumes (see next section). 



D. Finite size scaling and the critical temperature 



By definition, the correlation length characterizing the non-local degree of correlation of a system diverges at a 
critical point. Moreover, close enough to the transition it diverges as 

fcorr « lA-" (IV.7) 

where t = 1 — T/Tc and is a universal critical exponent. For the U{1) universality class (which contains superfluid 
phase transitions), this exponent is well-known: u = 0.671. 

When dealing with systems that have a finite size L^, the theory of the renormalization group (RG) predicts a 
very specific behavior for the correlation functions close enough to the transition temperature (see e.g. Ref. |38]). 
In particular, the two-body density matrix K{L, T) that gives the order parameter for off-diagonal long-range order, 
scales as 

R{L, T) = L^+'^KiL, T) ^ f{x){l + cL-'^ + ...) (IV.8) 

where 77 — 0.038 is another universal critical exponent, f{x) is a universal analytic function, x = {L / S^corrY ^ , and c 
is a non-universal constant, and lu ~ 0.8 is the critical exponent of the leading irrelevant field. One should keep in 
mind that typically one knows neither c nor Tc, but is interested in finding the latter. 

In a typical Monte Carlo calculation K{L, T) is computed for various lengths L,; and temperatures T. The procedure 
to locate the critical point (characterized by scale invariance) involves finding the "crossing" temperatures T!^ , for 
which R{Li,Tij) = R{Lj,Tij) at two given lengths Li and Lj. Assuming that one is close to the transition (so that 
the correlation length is large compared to any other scale), one can expand f{x{\t\)) — /(O) + f {0)L^^'^b\t\ (where 
^corr — b\t\~'^ was used), and derive the relation 

\T,-T,^\ = Kgm,Lj) (IV.9) 

where 



1 - 



&) -1 

1/^ 



(IV. 10) 



and K = cTcf{0)/bf'{0). If there were no non-universal corrections to scaling (i.e. if c = 0), then k = and Tc — Tij, 
which means that, upon scaling by the appropriate factor (as above) all the curves K{L, T) corresponding to different 
L's would cross exactly at Tc- In general these corrections are present, and it is therefore necessary to perform a 
linear fit of Tij vs. g{Li, Lj) and extrapolate to infinite L in order to determine the true Tc. Following such procedure 
our data for the condensate fraction of the unitary Fermi gas indicates that Tc ^ 0.15(l)eF, considerably lower than 
the characteristic temperature Tg = 0.23(2) found by studying the behavior of the energy and the chemical potential. 
Even though this result for Tc is close to estimates by other groups (see e.g. [3|)i it should be pointed out that the 
experimental data of Ref. [s^l shows a distinctive feature in the energy versus entropy curve at a temperature close 
to To (see Ref. {35]). whereas a clear signature of a transition at a lower temperature remains to be found. 



V. RESULTS AWAY FROM UNITARITY 



In the following we describe the results of our calculations away from unitarity. The system was placed on a lattice of 
volume V — (8/)'^, filled with = 45± 15 particles. In all cases the temperatures cover the range 0.07 < T/ep < 0.5, 
corresponding to Nr steps in the imaginary time direction varying from A^,- ~ 1700 to A^t- ~ 200, respectively. 
The temperature is limited from below by the precision of our computers (because the matrices involved become 
ill-conditioned in the sense explained in section III B), and from above by the fact that our phase space has a natural 
UV cutoff given by the inverse lattice spacing. In all cases the occupation of the high-energy modes smaller than 
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FIG. 11: (Color online) Condensate fraction a{T) scaled with the lattice size L, black squares for 6'^, red stars for 8'' and blue 
circles for 10^ (highest on the left) lattices, respectively for the case l/kpa — —0.1. The errorbars correspond to the statistical 
errors. 

1% percent. The coupling strength was varied in the range —0.5 < 1/kpa < 0.2 (where kp = (37r^n)^/^), and is 
limited on the negative (BCS) side by the finite volume V (which may become comparable to the size of the Cooper 
pairs deep in the BCS regime), and on the positive (BEC) side by the finite lattice spacing I (whose size eventually 
becomes inadequate to describe localized dimers of size a — 0{l), deep in the BEC regime, and which manifests 
itself as poor convergence of observables). The number of uncorrelated Monte Carlo samples varies from 7500 at the 
lowest temperatures to 2500 at the highest. The Monte Carlo auto-correlation time was ~ 200 samples (estimated by 
studying the autocorrelation of the energy), implying a statistical error of less than 2%. 
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0.24(8) 
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0.56(8) 


0.17(1) 


0.35(1) 
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0.26(3) 
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TABLE I: Results for the ground state energy, the characteristic temperature To, and the corresponding chemical potential 
and energy, from the caloric curves of Fig. 1121 and the upper bounds on the critical temperature Tc from finite size scaling and 
the corresponding chemical potential and energy. 

Repeating the analysis of Sec lIVDl one arrives at the estimates for Tc shown in Table HI Notice however that, given 
the rather small lattice sizes used (the limitation being given by the required computer power/time, and ultimately 
by its scaling with the size of the lattice), an extrapolation to L ^ oo is difficult. Still, the study of the crossing 
temperatures provides us at least with upper bounds on which is what we show in Table HI Unfortunately, 
it was not possible for us to explore temperatures below O.le^i and therefore we were unable to find T^ on the BCS 
side beyond [kpo)^^ = —0.3. The fact that our upper bound on Tc at unitarity agrees with extrapolations to T — > oo 
performed by other groups 0] indicates that these bounds are not far from the actual result. In Ref. the authors 
performed a finite size scale analysis of our initial data [B| and found a value of Tc in agreement with their result. In 
Fig. [Tl]we show an example of the finite size scaling analysis at 1/kpa = —0.1 performed as explained in Sec lIVDI 
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FIG. 12; (Color online) Energy (left panel) and chemical potential (right panel) as functions of T jep and l/kpa. The dashed 
line shows the location of Tc and the solid line represents To, both as functions of l/kpa, from Table H] 

The last two columns of Table U show the chemical potential and the energy at the value of the bound on the critical 
temperature. In particular, the values at unitarity, namely {kpa)~^ — Q;^c/£f — 0-4:3{l); Ec/ Ep = 0.45(1) should be 
compared with the results of reference 0: fJ'c/sp = 0.493(14); i?c/-E_F — 0.52(1), both of which are higher than our 
estimates. It should also be pointed out that the latter values, shown in Fig. [151 slightly violate the bounds imposed 
by thermodynamic stability (see Appendix), which is not the case for our data. 

Recently the Amherst-ETH group has posted values Tc for a couple of values of the coupling constant l/kpa > 
with some of the details of the calculations, see Ref. [l^ and Fig. [O] 

The strong dependence observed by these authors earlier on the filling factor, see discussion around Fig. ([3]), was 
presumably due to the use of an inaccurate representation of the kinetic energy, which becomes accurate only at very 
low fiUing factors (when kp < 1/0; see, however, our discussion on density dependence in Sec. IIIIDI The values 
of the critical temperature estimated in this work and in Ref. [221 agree within the error bars at unitarity and at 
1/kpa « 0.2. The value Tc/ep = 0.252(15) at 1/kpa = 0.474(8) [22| does not seem to follow the systcmatics suggested 
by the rest of the results for 1/kpa < 0.22, the critical temperature for the hard and soft bosons ^SQ] and the limiting 
BCS and BEC behavior. If one ignores the value Tc/ep = 0.252(15) at 1/kpa — 0.474(8) [111, the data presented in 
Fig. [T3[ would thus suggest that for the value of the coupling constant 1/kpa « 0.8 the critical temperature attains a 
maximum of Tc/ep « 0.23(2). 

The results for the energy E per particle (in units of the free gas ground-state energy Ep — 3/5epN, where N 
is the total number of particles) are shown in Fig. I12[ along with the chemical potential /i (in units of the free gas 
Fermi energy ep). For every value of 1/kpa that we studied, our data presents two salient features: below certain 
temperature Tq, fi/ep is approximately constant (a feature of free Bose gases in the condensed phase); above that 
temperature fj,/ep decreases steadily, while E/Ep becomes the energy of a free Fermi gas, offset by a constant energy 
(whose specific value depends on 1/kpa). In Table[Tlwe summarize our results at Tq. The errors represent uncertainties 
in the point of departure of E/Ep from the (offset) free Fermi gas, and the departure of fj,/ep from its (approximately) 
constant low-temperature value. In the same table we also show our extrapolated values for the ground state energies. 
At low enough temperature both E/Ep and fi/ep become approximately constant. From this observation it can 
be inferred that the thermal fluctuations are small enough that those constant values should not differ greatly from 
their ground-state values. For reference we also include our data for the system at unitarity. The fact that the latter 
data falls in the right place shows that our calculations are quite close to the dilute limit. In this respect one should 
also note that the T = fixed-node Monte-Carlo calculations show a similar agreement, even though in one case 
nrj) « 10"'^ 0, while in the other nrQ « 10~^ where is the effective range of the interaction used. Notice that 
our estimated value for ^ at unitarity is lower than the variational estimates reported in Refs. [3, H, [3, [l3| and in 
apparent agreement within error bars with the unpublished results of Ref. [2p| . 

In the left panel of Fig. [TH we show E/Ep extrapolated to T = 0, as a function of 1/kpa, together with the 
ground-state energy as determined in Refs. 0, [^, [l3|- Since our calculations are in principle exact, it is not suprising 
that our extrapolations yield results that are consistently lower than those in Refs. [1, |3| , since those calculations 
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BCS/BEC limits 



-?.5 -1 -0.5 0.5 1 1.5 

FIG. 13: (Color online) The solid (purple) curve and diamonds represent the temperature To at which the energetic behavior 
of the system shows a transition from Bose-like to Fermi-like behavior. The solid (green) curve in the low-left corner is Tc 
estimated within the BCS theory with the Gorkov-Melik-Barkhudarov correction [2^. The solid (green) curve on the right 
is the Tc « 0.218eF for a non-interacting Rose gas. The (red) up-triangles and the down-triangles joined by a dotted line 
respectively are the Tc for the hard-core and soft-core bosons calculated in Ref. [s^. The three (blue) dots with error bars 
joined by a dashed line are the results of Ref. ^2^], while the six (blue) squares joined by a solid line are our new estimates for 
Tc presented din this work. 



are based on the Fixed-Node approximation, which are thus variational and provide an estimate for the ground state 
energy from above. 



ergy irom above. 

In the right panel of Fig. [14] we present the results for the pairing gap at the lowest temperatures. This quantity 
IS determined via a calculation of a response function x (see Ref. |40| ) as a function of momentum p: 



was 



where 



X(p)-- / rfr^Mp.r) (V.l) 



Tr[e-(^'-^)(^-^^) (p)e-^(-f^-^^) V| (p)] 



^Mp, r) = ' ^ (V.2) 

z(/3,/i,y) 

is the temperature Green function. This response function has been shown by us in Ref. [40l | to be accurately 
parametrized by the independent quasi-particle form given by 

1 e^^^^P) - 1 

with 



i.(p) = ^(g + C/-,)%A^. (V.4) 

In this expression, a = m/m,,^ where is an effective mass, U is the mean-field potential, A represents the pairing 
gap, and ^ is the chemical potential. All of these quantities are functions of temperature (see Ref. [i^ for further 
details). The data for A shown in the right panel of Fig. [T4l corresponds to the lowest temperatures we have simulated 
(namely T < O.lep). Our results for A agree qualitatively with the data by other groups ([1, [l^, [ij)- Away from 
unitarity, however, our data falls systematically below the data by other groups. 
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FIG. 14: (Color online) Ground state energy (left panel, in units of the free gas ground state energy S/SefN), obtained by 
extrapolating our finite temperature results to T = 0, is shown in red circles with error bars (representing the uncertainty in 
the extrapolation). The data of Ref. 3j appears in blue circles, Ref. in black dash-dotted line, and Ref. 4] in green squares. 
The right panel shows the gap at the lowest temperatures, from our work in Ref. jlQ], in red circles. Blue circles with error 
bars show the data of Ref. Q, the black dash-dotted line shows the data of Ref. and the green triangle represents the data 
of Ref 



VI. SUMMARY AND CONCLUSIONS 



In this paper we have described the technical details involved in the non-perturbative calculation of thermal averages 
of systems of interacting fermions at finite temperature. We have performed calculations of the thermal properties of 
a system of spin 1/2 fermions at and away from the unitary point. The particles were placed on a 3D spatial lattice, 
in a path integral formulation of the interacting many-body problem. 

By studying the finite size scaling of the condensate fraction we have established upper bounds on the critical 
temperature of the superfluid-normal phase transition, for couplings around the unitary point in the region —0.5 < 
(/cfi)"^ < 0.2. At unitarity we find Tc ~ 0.15(1), which is in agreement with Ref. m. In contrast, at the transition 
the energy Ec and the chemical potential fic that we find are lower than those of Ref. [31 by about 15%. Furthermore, 
we find that Ec and fic of Ref. [Tj] slightly violate the bounds imposed by thermodynamic stability, which are satisfied 
by our data, as shown in Fig. [151 

For all the couplings we studied, in particular at unitarity, our results for the universal function ^ and the chemical 
potential are consistent with normal Fermi gas behavior above a characteristic temperature To > Tc that depends on 
the coupling. Tq is obtained by studying the deviations of the caloric curve from that of a free Fermi gas. Furthermore, 
the chemical potential is approximately constant below Tq. The existence of such a characteristic temperature that 
is different from the critical temperature is analogous to the case of water, where density reaches a maximum at a 
temperature T ~ 4°C, which is above the T = 0°C liquid-solid phase transition. 

At unitarity we find Tq = 0.23(2), which is in agreement with the experimental results of Ref. [s^l, where measure- 
ments of the caloric curve and energy vs. entropy curve of a unitary Fermi gas were reported. For Tc < T < Tq there 
is a noticeable departure from normal Fermi gas behavior, possibly due to pairing effects. 

Extrapolations of our data for the energy to T = are systematically below the results by other groups. This is 
not surprising because Green Function Monte Carlo methods provide an upper bound to the energy. 

We also compare our low temperature results for the gap (determined through the calculation of a response function, 
as explained in Ref. j4(?|) with ground state calculations and find reasonably good agreement close to the unitary point, 
and somewhat lower values on the BCS side of the resonance. 



VII. APPENDIX: THERMODYNAMIC RELATIONS AT UNITARITY. 



In this section we complete our discussion of thermodynamics at unitarity by deriving a number of identities and 
expressing the various thermodynamic functions in useful forms. 

We start with the grand-canonical ensemble, where the thermodynamics is derived from the thermodynamic poten- 
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tial 0,{T, n, V) = —PV. At unitarity, where both n and T are conventionally measured in units of the free gas Fermi 
energy ep-, we can write in terms of a function hx{z), where z = fi/T, as in Ref. [35| : 

n{T, V) = n{z, V) = -P{z)V - -1/3 [ThTiz)f^ V (VII.l) 

where /? — (7^^)'^^^- This form is useful because thermodynamic stability implies three conditions on Ht'- hx > 0, 
h'rp > Q and hi^ > 0. The form of this function is shown in Fig. [TSl along with data of Ref. [TJ- 




FIG. 15: (Color online) MC data from (blue circles), Ref. 0] (six black points). The four straight lines starting at the origin 
are the T ^ limits of hriz 00) = 2^''^2/^?''^ where = 0.42(2) 0, IH, = 0.59 for meanfield/BCS approximation 
and 5 = 1 for the free Fermi gas model respectively. The two solid lines (red/lower and green/higher) correspond to /it(z) 
calculated in the free Fermi gas and the BCS/meanfield approximation hriz) respectively. 

The particle number density and the entropy per particle can then be derived as follows: 



N 
V 



1 dfl 



5Ph'j.{z) 
2ThT{z) 



(VIL2) 



S_ 

N 



1 dVL 
Ndf 



5_P_ 
2^ 



h'j.{z) 



hriz) 



hT(z) 



h'j,{z) 



(VII.3) 



Using the thermodynamic identity E = TS — PV + fiN , inserting the expressions above, we find 



E 



-PV 



h'j,{z) 



hriz) 



5 hUz) 



-PV - PV 



-PV 



(VII.4) 



Using this relation together with the constant volume identity dE/dT — TdS/dT, one can derive relation pil.9p . 
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